Comment on ‘ ‘ On using a too large integration time step in molecular dynamics simulations of coarse - grained molecular models

In a recent study published in this journal, Winger et al. investigate the magnitude of the time step to integrate the equations of motion in simulations with the coarse-grained MARTINI force field, using an implementation of MARTINI in the GROMOS software. Based primarily on the drift in temperature and the magnitude of the energy fluctuations in bulk liquids, the authors conclude that a time step not exceeding 10 fs should be used to avoid artificial energy flow into or out of the system. In most applications of the MARTINI model to date, time steps of 20–40 fs have been used. Thus, the observation of Winger et al. raises questions about possible artefacts caused by the apparent use of a too large integration time step. Although we appreciate the effort put into testing our model, we do not support the conclusions drawn. First, the MARTINI force field is not an atomistically detailed force field. Many assumptions underlie the model, the major one being the neglect of some of the atomistic degrees of freedom. As a result, the interactions between particles are effective ones and the energy landscape is highly simplified. This simplified energy landscape allows for a greatly increased sampling speed at the cost of a loss of detail. This makes CG models in general so powerful. The emphasis, therefore, should not be to sample the energy landscape as accurately as possible, but rather, as effectively as possible. This is in contrast to traditional all-atom models, for which the energy landscape is more realistic and an accurate integration scheme is more important. In practice, the inherent ‘fuzziness’ of the MARTINI model makes the presence of artificial small energy sinks or sources a less critical problem than in accurate atomistic simulations. Most importantly, structural properties are very robust with respect to time step; Winger et al. show that, even for the worst-case scenario (50 fs time step), there are no noticeable effects on structural properties of the systems investigated. Moreover, thermodynamic properties such as the free energy of solvation also appear insensitive to the size of the time step. Thus, if the goal is to generate representative ensembles quickly, large time steps seem acceptable. Secondly, using the standard GROMACS implementation of the MARTINI model, we are unable to reproduce the results of Winger et al. Analogous to the systems investigated by Winger et al., we simulated bulk water and hexadecane systems, both at NVE and NpT conditions. Our systems are composed of 1700 CG water particles or 800 hexadecane molecules. In the NVE set-up, the box volume was 205.5 nm (water) and 362 nm (hexadecane). For the NpT ensemble, weak coupling to a temperature bath (T = 300 K, coupling constant tT = 1 ps) and pressure bath (p = 1 bar, tp = 5 ps, compressibility 5 10 4 bar ) was used. The standard MARTINI protocol for the treatment of non-bonded interactions was followed, i.e. using the shifted potential F as implemented in GROMACS: